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The Phase Transition of the Spin-1/2 Heisenberg Model with a Spatially Staggered 

Anisotropy on the Square Lattice. 

F.-J. Jiang^El 

1 Center for Theoretical Physics, Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139 

Puzzled by the indication of a new critical theory for the spin-1/2 Heisenberg model with a 
spatially staggered anisotropy on the square lattice as suggested in [J], we re-investigate the phase 
transition of this model induced by dimerization. We focus on studying the finite-size scaling of 
the observables p s \L and PsiT, where L stands for the spatial box sizes used in the simulations 
and psi with i £ {1, 2} is the spin-stiffness in i-direction. We find by performing finite-size scaling 
using the observable p s iL, which corresponds to the spatial direction with a fixed antiferromagnetic 
coupling, one would suffer a much less severe correction compared to that of using p s \L. Therefore 
psiL is a better quantity than p s \L for finite-size scaling analysis concerning the limitation for the 
availability of large volumes data in our study. Remarkably, by employing the method of fixing the 
aspect-ratio of spatial winding numbers squared in the simulations, even from p B \L which receives 
the most serious correction among the observables considered in this study, we arrive at a value for 
the critical exponent v which is consistent with the expected 0(3) value by using only up to L — 64 
data points. 

PACS numbers: 



I. INTRODUCTION 



Heiscnberg-type models have been studied in great de- 
tail during the last twenty years because of their phe- 
nomenological importance. For example, it is believed 
that the spin-1/2 Heisenberg model on the square lat- 
tice is the correct model for understanding the undoped 
precursors of high T c cuprates (undoped antiferromag- 
nets). Further, due to the availability of efficient Monte 
Carlo algorithms as well as the increasing power of com- 
puting resources, properties of undoped antiferromagnets 
on geometrically non-frustrated lattices have been deter- 
mined to unprecedented accuracy 0, H, 0, [H, @, 0, Hi- 
For instance, using a loop algorithm, the low-energy pa- 
rameters of the spin-1/2 Heisenberg model on the square 
lattice are calculated very precisely and are in quanti- 
tative agreement with the experimental results De- 
spite being well studied, several recent numerical investi- 
gation of anisotropic Heisenberg models have led to un- 
expected results [l], [Io|, EH- I n particular, Monte Carlo 
evidence indicates that the anisotropic Heisenberg model 
with staggered arrangement of the antiferromagnetic cou- 
plings may belong to a new universality class, in contra- 
diction to the theoretical 0(3) universality prediction [l[. 
For example, while the most accurate Monte Carlo value 
for the critical exponent v in the 0(3) universality class 
is given by v = 0.7112(5) |l2|, the corresponding v de- 
termined in [l[ is shown to be v = 0.689(5). Although 
subtlety of calculating the critical exponent v from per- 
forming finite-size scaling analysis is demonstrated for a 
similar anisotropic Heisenberg model on the honeycomb 



lattice [13j, the discrepancy 
v = 0.7112(5) observed in 
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between v — 0.689(5) and 
remains to be understood. 
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In order to clarify this issue further, we have simulated 
the spin-1/2 Heisenberg model with a spatially staggered 
anisotropy on the square lattice. Further, we choose to 
analyze the finite-size scaling of the observables p s iL and 
p S 2L, where L refers to the box sizes used in the simu- 
lations and p s i with i g {1,2} is the spin stiffness in 
i-direction. The reason for choosing p s \L and p S 2L is 
twofold. First of all, these two observables can be cal- 
culated to a very high accuracy using loop algorithms. 
Secondly, one can measure p s \ and p s i separately. In 
practice, one would naturally either measure p s which 
is the average of p s \ and p S 2 in order to increase the 
statistics, or p S 2, which corresponds to the spatial direc- 
tion with a fixed antiferromagnetic coupling, would not 
be used for data analysis since more measurements is re- 
quired in order to obtain a good statistics for this observ- 
able. However for the model considered here, it is useful 
to measure quantities which are sensitive to anisotropy. 
Surprisingly, as we will show later, the observable p S 2L 
receives a much less severe correction than p s \L does. 
Hence p S 2L is a better observable than p s \L (or p s L) for 
finite-size scaling analysis concerning the limitation for 
the availability of large volumes data in this study. In 
addition, instead of using a fixed aspect-ratio of spatial 
box sizes as done in most Monte Carlo calculations, in our 
investigation we employ the method of fixing the aspect- 
ratio of spatial winding numbers squared in the simula- 
tions which we will introduce briefly later. Remarkably, 
combining the idea of fixing the aspect-ratio of spatial 
winding numbers squared in the simulations and finite- 
size scaling analysis, unlike the unconventional value for 
v observed in [l[, even from p s \L which suffers a very 
serious correction, we arrive at a value for v which is 
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FIG. 1: The anisotropic Heisenberg model considered in this 
study. 

consistent with that of 0(3) by using only up to L = 64 
data points. 

This paper is organized as follows. In section [Til the 
anisotropic Heisenberg model and the relevant observ- 
ables studied in this work are briefly described. Section 
IIIII contains our numerical results. In particular, the cor- 
responding critical point as well as the critical exponent 
v are determined by fitting the numerical data to their 
predicted critical behavior near the transition. Finally, 
we conclude our study in section ITVl 

II. MICROSCOPIC MODEL AND 
CORRESPONDING OBSERVABLES 

The Heisenberg model considered in this study is de- 
fined by the Hamilton operator 

H=Y, JS x -S y + J 'Sx' -Sy>, (1) 

(xy) {x'y') 

where J' and J are antiferromagnetic exchange couplings 
connecting nearest neighbor spins (xy) and (x'y 1 ), re- 
spectively. Figure 1 illustrates the Heisenberg model de- 
scribed by Eq. ([JJ). To study the critical behavior of this 
anisotropic Heisenberg model near the transition driven 
by the anisotropy, in particular to determine the critical 
point as well as the critical exponent v, the spin stiff- 
nesses in the 1- and 2-directions which are defined by 

Psi = ^<W?>. (2) 

are measured in our simulations. Here (3 is inverse tem- 
perature and L refers to the spatial box sizes. Further 
W? is the winding number squared in the z-direction. By 
carefully investigating the spatial volumes and the J'/J 
dependence of p S iL, one can determine the critical point 
as well as the critical exponent v with high precision. 

III. DETERMINATION OF THE CRITICAL 
POINT AND THE CRITICAL EXPONENT v 

To calculate the relevant critical exponent v and to 
determine the location of the critical point in the pa- 
rameter space J' I J, one useful technique is to study the 
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FIG. 2: p s L (upper panel) and p a iL (lower panel) as functions 
of J' I J. 

finite-size scaling of certain observables. For example, if 
the transition is second order, then near the transition, 
the observable p sl L for i G {1,2} should be described 
well by the following finite-size scaling ansatz 

L {t) = {l+bL-")go{tL 1 '»), (3) 

where Ol stands for p S iL, t = (j c —j)/je with j — (J'/ J), 
b is some constant, v is the critical exponent correspond- 
ing to the correlation length £ and u is the confluent 
correction exponent. Finally go appearing above is a 
smooth function of the variable tL 1 ^. From Eq. ©, one 
concludes that the curves of different L for Ol, as func- 
tions of J' j J, should have the tendency to intersect at 
critical point (J'/J) c for large L. In the following, we 
will employ the finite-size scaling formula, Eq. Q, for 
PsiL with i G {1,2} to calculate the critical exponent v 
and the critical point (J'/J) c . Without losing the gen- 
erality, in our simulations we have fixed J to be 1.0 and 
have varied J' . Further, the box size used in the simu- 
lations ranges from L — 6 to L = 64. We also use large 
enough j3 so that the observables studied here take their 
zero-temperature values. Figure [2] shows the observables 
p s L and p S 2L as functions of J' /J. The figure clearly 
indicates the phase transition is second order since dif- 
ferent L curves for both p s L and p s -2,L tend to intersect 
at a particular point in the parameter space J' / J. What 
is the most striking observation from our results is that 
the observable p s L receives a much severe correction than 
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FIG. 3: Fits of p s L (upper panel) and p s iL (lower panel) 
to Eq. [3] While the circles and squares on these two panels 
are the numerical Monte Carlo data from the simulations, the 
solid curves are obtained by using the results from the fits. 



p s %L docs. This can be understood from the trend of the 
crossing among these curves of different L in figure [2] 
Therefore one expects a better determination of v can be 
obtained by applying finite-size scaling analysis to p S 2L. 
Before presenting our results, we would like to point out 
that since data from large volumes might be essential 
in order to determine the critical exponent v accurately 
as suggested in [l3| . we will use the strategy employed 
in [13[ for our data analysis as well. A Taylor expan- 
sion of Eq. ([3]) up to fourth order in tl}l v is used to fit 
the data of p S 2L. The critical exponent v and critical 
point (J'/J) c calculated from the fit using all available 
data of p a iL are given by 0.6934(13) and 2.51962(3), re- 
spectively. The upper panel of figure [3] demonstrates the 
result of the fit. Notice both v and ( J'/ J) c we obtain are 
consistent with the corresponding results found in [l[ . By 
eliminating some data points of small L. we can reach a 
value of 0.700(3) for v by fitting p s ^L with L > 26 to 
Eq. ((3]). On the other hand, with the same range of L 
(L > 26), a fit of p s L to Eq. [3] leads to v = 0.688(2) 
and (J'/J) c = 2.5193(2), both of which are consistent 
with those obtained in [1[ as well (lower panel in figure 
[3]). By eliminating more data points of p s L with small 
L, the values for v and (J'/J) c calculated from the fits 



are always consistent with those quoted above. What 
we have shown clearly indicates that one would suffer 
the least correction by considering the finite-size scaling 
of the observable p S 2L. As a result, it is likely one can 
reach a value for v consistent with its 0(3) prediction, 
namely v = 0.7112(5) if large volume data points for p S 2 
are available. Here we do not attempt to carry out such 
task of obtaining data for L > 64. Instead, wc employ 
the technique of fixing the aspect-ratio of spatial winding- 
numbers squared in the simulations. Surprisingly, com- 
bining the idea of fixing the aspect-ratio of winding num- 
bers squared and finite-size scaling analysis, even from 
the observable p a \L which is found to receive the most 
severe correction among the observables considered here, 
we reach a value for the critical exponent v consistent 
with v = 0.7112(5) without additionally obtaining data 
points for L > 64. The motivation behind the idea of fix- 
ing the aspect-ratio of spatial winding numbers squared 
in the simulations is as follows. Intuitively the winding 
numbers squared W 2 and W% indicate the ability of the 
loops moving around 1- and 2-directions, respectively. 
Further, one can consider the original anisotropic system 
on the square lattice as an isotropic system on a rectangu- 
lar lattice. From these points of view, it is (W 2 )/(yV%) : 
not (L2/L1) 2 , plays the role of the quantity (L^/L\) 2 
for the system, here Li and L\ with i G {1,2} are the 
spatial box size used in the simulations and the linear 
physical length of the system in i-direction, respectively. 
Indeed it is demonstrated in Q that rectangular lattice 
is more suitable than square lattice for studying the spa- 
tially anisotropic Hcisenbcrg model with different antifer- 
romagnctic couplings J\, J2 in 1- and 2-directions. The 
idea of fixing the aspect-ratio of spatial winding numbers 
squared quantifies the method used in 0]. In general for 
a fixed Z/2, one can vary L\ and J'/J in order to reach 
the criterion of a fixed aspect-ratio of spatial winding 
numbers squared in the simulations. For our study here, 
without obtaining additional data, this method is imple- 
mented as follows. First of all, we calculate (W 2 ) / (Wf) 
for the data point at J' / J = 2.5196 with L = 40 which we 
denote by w/. Notice since only the aspect-ratio of the 
linear physical lengths squared is fixed, we choose L\ = L 
for our finite-size scaling analysis. After obtaining this 
number, a linear interpolation for other data points of p s \ 
based on (w/wf)^ 1 ^ 2 ^ is performed in order to reach the 
criterion of a fixed aspect-ratio of spatial winding num- 
bers squared in the simulations. The w appearing above 
is the corresponding (W 2 )/{W 2 ) for other data points. 
Further, we keep the number \w/wf — 1| smaller than 
0.055 so that the interpolation results are reliable. A fit 
of the interpolated (p s \)- m L data to Eq. [3] with lu being 
fixed to its 0(3) value (lu = 0.78) leads to v = 0.706(7) 
and {J/J) c = 2.5196(1) for 36 < L < 64 (figure [4}. The 
subscript " in " appearing above stands for interpola- 
tion. Letting to be a fit parameter results in consistent 
v = 0.707(8) and (J'/J) c = 2.5196(7). Further, we al- 
ways arrive at consistent results with v = 0.706(7) and 
(J'/J) c = 2.5196(1) from the fits using L > 36 data. The 
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(J7J) = 2.5196(1) 
v = 0.706(7) 
X 2 /d.o.f. = 0.5 




FIG. 4: Fit of interpolated (p sl ) in L data to Eq.0 While the 
circles are the numerical Monte Carlo data from the simula- 
tions, the solid curves are obtained by using the results from 
the fit. 



value of v we calculate from the fit is in good agreement 
with the expected 0(3) value v = 0.7112(5). The critical 
point (J'/J) c = 2.5196(1) is consistent with that found 
in [l[ as well. To avoid any bias, we perform another 
analysis for the same set of Monte Carlo data without 
interpolation. By fitting this set of original data points 
to Eq. [3] with a fixed uj = 0.78, we arrive at v = 0.688(7) 
and (J'/J) c = 2.5197(1) (figure ©, both of which again 
agree quantitatively with those determined in pj. Fi- 
nally we would like to make a comment regarding the 
choice oi Wf. In principle one can calculate wj for any 
L and for any J'/J close to the critical point. How- 
ever since it is shown in [l3| that data points of larger 
volumes is essential for a quick convergence of v, it will 
be desirable to choose Wf such that the set of interpo- 
lated data contains sufficiently many data points from 
large volumes. For example, using the Wf obtained at 
J'/J = 2.5191 (J'/J = 2.5196) with L = 40 (L = 44), we 
reach the results of v = 0.704(7) and (J'/J) c = 2.5196(1) 
(u = 0.705(7) and (J'/J) c = 2.5196(1)) from the fit 
with a fixed u) = 0.78. These values for v and (J'/J) c 
agree with what we have obtained earlier. Interestingly, 
it seems that the idea of fixing the aspect-ratio of wind- 
ing numbers squared determines the critical exponent v 



more accurately than the conventional method of fixing 
the aspect-ratio of box sizes in the simulations. It would 
be interesting to explore this new method systematically 
including applying it to the study of phase transition for 
other spatially anisotropic Heisenberg models. 



IV. DISCUSSION AND CONCLUSION 

In this letter, we revisit the phase transition of the 
spin-1/2 Heisenberg model with a spatially staggered 
anisotropy. We find that the observable ps^L suffers a 
much less severe correction compared to that of p s \L, 
hence is a better quantity for finite-size scaling analysis. 
Further, by employing the method of fixing the aspect- 
ratio of spatial winding numbers squared in the simula- 
tions, we arrive at v = 0.706(7) for the critical exponent v 
which is consistent with the most accurate Monte Carlo 
0(3) result v = 0.7112(5) by using only up to L = 64 
data points derived from p s \L. 

The simulations in this study were performed using 
the ALPS library 14| on personal desktops. This work 
is supported in part by funds provided by the DOE Office 
of Nuclear Physics under grant DE-FG02-94ER40818. 
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FIG. 5: Fit of original p s iL data to Eq. [3] While the circles 
are the numerical Monte Carlo data from the simulations, the 
solid curves are obtained by using the results from the fit. 
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